Library of Basic Chemical Reaction Network Models

Below we will present various simple and established chemical reaction network (CRN) models. Each model is given some brief background, implemented using the @reaction_network DSL, and basic simulations are performed.

Birth-death process

The birth-death process is one of the simplest possible CRN models. It consists of a single component ($X$) which is both produced and degraded at linear rates:

using Catalyst
bd_process = @reaction_network begin
    (p,d), ∅ <--> X
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]

Next, we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).

u0 = [:X => 1]
tspan = (0.0, 10.0)
ps = [:p => 1, :d => 0.2]

We can now simulate our model using all three interpretations. First, we perform a reaction rate equation-based ODE simulation:

using OrdinaryDiffEq
oprob = ODEProblem(bd_process, u0, tspan, ps)
osol = solve(oprob, Tsit5())

Next, a chemical Langevin equation-based SDE simulation:

using StochasticDiffEq
sprob = SDEProblem(bd_process, u0, tspan, ps)
ssol = solve(sprob, ImplicitEM())

Next, a stochastic chemical kinetics-based jump simulation:

using JumpProcesses
dprob = DiscreteProblem(bd_process, u0, tspan, ps)
jprob = JumpProblem(bd_process, dprob, Direct())
jsol = solve(jprob, SSAStepper())

Finally, we plot the results:

using Plots
oplt = plot(osol; title = "Reaction rate equation (ODE)")
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
plot(oplt, splt, jplt; lw = 3, size=(800,700), layout = (3,1))

Michaelis-Menten enzyme kinetics

Michaelis-Menten enzyme kinetics is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions, it can be simplified to a single function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:

using Catalyst
mm_system = @reaction_network begin
  kB, S + E --> SE
  kD, SE --> S + E
  kP, SE --> P + E
end

\[ \begin{align*} \mathrm{S} + \mathrm{E} &\xrightleftharpoons[kD]{kB} \mathrm{SE} \\ \mathrm{SE} &\xrightarrow{kP} \mathrm{P} + \mathrm{E} \end{align*} \]

Next, we perform ODE, SDE, and jump simulations of the model:

u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
tspan = (0., 100.)
ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]

using OrdinaryDiffEq
oprob = ODEProblem(mm_system, u0, tspan, ps)
osol  = solve(oprob, Tsit5())

using StochasticDiffEq
sprob = SDEProblem(mm_system, u0, tspan, ps)
ssol = solve(sprob, ImplicitEM())

using JumpProcesses
dprob = DiscreteProblem(mm_system, u0, tspan, ps)
jprob = JumpProblem(mm_system, dprob, Direct())
jsol = solve(jprob, SSAStepper())

using Plots
oplt = plot(osol; title = "Reaction rate equation (ODE)")
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1))

SIR infection model

The SIR model is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery.

using Catalyst
sir_model = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
end

\[ \begin{align*} \mathrm{S} + \mathrm{I} &\xrightarrow{\alpha} 2 \mathrm{I} \\ \mathrm{I} &\xrightarrow{\beta} \mathrm{R} \end{align*} \]

First, we perform a deterministic ODE simulation:

using OrdinaryDiffEq, Plots
u0 = [:S => 99, :I => 1, :R => 0]
tspan = (0.0, 500.0)
ps = [:α => 0.001, :β => 0.01]

# Solve ODEs.
oprob = ODEProblem(sir_model, u0, tspan, ps)
osol = solve(oprob, Tsit5())
plot(osol; title = "Reaction rate equation (ODE)", size=(800,350))

Next, we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of an outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.

using JumpProcesses
dprob = DiscreteProblem(sir_model, u0, tspan, ps)
jprob = JumpProblem(sir_model, dprob, Direct())

jsol1 = solve(jprob, SSAStepper())
jsol2 = solve(jprob, SSAStepper())
jsol3 = solve(jprob, SSAStepper())

jplt1 = plot(jsol1; title = "Outbreak")
jplt2 = plot(jsol2; title = "Outbreak")
jplt3 = plot(jsol3; title = "No outbreak")
plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))

The Wilhelm model

The Wilhelm model was introduced in Wilhelm (2009) as the smallest CRN model (with constant rates) that exhibits bistability.

using Catalyst
wilhelm_model = @reaction_network begin
    k1, Y --> 2X
    k2, 2X --> X + Y
    k3, X + Y --> Y
    k4, X --> 0
end

\[ \begin{align*} \mathrm{Y} &\xrightarrow{k1} 2 \mathrm{X} \\ 2 \mathrm{X} &\xrightarrow{k2} \mathrm{X} + \mathrm{Y} \\ \mathrm{X} + \mathrm{Y} &\xrightarrow{k3} \mathrm{Y} \\ \mathrm{X} &\xrightarrow{k4} \varnothing \end{align*} \]

We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states.

using OrdinaryDiffEq, Plots
u0_1 = [:X => 1.5, :Y => 0.5]
u0_2 = [:X => 2.5, :Y => 0.5]
tspan = (0., 10.)
ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]

oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
osol1 = solve(oprob1, Tsit5())
osol2 = solve(oprob2, Tsit5())
plot(osol1; lw = 4, idxs = :X, label = "X(0) = 1.5")
plot!(osol2; lw = 4, idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,350))

Simple self-activation loop

The simplest self-activation loop consists of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.

using Catalyst
sa_loop = @reaction_network begin
    v₀ + hill(X,v,K,n), ∅ --> X
    d, X --> ∅
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{v_0 + \frac{v X^{n}}{K^{n} + X^{n}}} \mathrm{X} \end{align*} \]

A simple example of such a loop is a transcription factor which activates its own gene. Here, $v₀$ represents a basic transcription rate (leakage) in the absence of the transcription factor.

We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.

using JumpProcesses, OrdinaryDiffEq, Plots
u0 = [:X => 4]
tspan = (0.0, 1000.0)
ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]

oprob = ODEProblem(sa_loop, u0, tspan, ps)
osol = solve(oprob, Tsit5())

dprob = DiscreteProblem(sa_loop, u0, tspan, ps)
jprob = JumpProblem(sa_loop, dprob, Direct())
jsol = solve(jprob, SSAStepper())

plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide="X", size = (800,350))

The Brusselator

The Brusselator is a well-known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").

using Catalyst
brusselator = @reaction_network begin
    A, ∅ --> X
    1, 2X + Y --> 3X
    B, X --> Y
    1, X --> ∅
end

\[ \begin{align*} \varnothing &\xrightarrow{A} \mathrm{X} \\ 2 \mathrm{X} + \mathrm{Y} &\xrightarrow{1} 3 \mathrm{X} \\ \mathrm{X} &\xrightarrow{B} \mathrm{Y} \\ \mathrm{X} &\xrightarrow{1} \varnothing \end{align*} \]

It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when combinatorial adjustment of rates is not performed. Since Catalyst automatically perform these adjustments, and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.

using OrdinaryDiffEq, Plots
u0 = [:X => 1.0, :Y => 1.0]
tspan = (0., 50.)
ps1 = [:A => 1.0, :B => 1.0]
ps2 = [:A => 1.0, :B => 1.8]

oprob1 = ODEProblem(brusselator, u0, tspan, ps1)
oprob2 = ODEProblem(brusselator, u0, tspan, ps2)
osol1  = solve(oprob1, Rodas5P())
osol2  = solve(oprob2, Rodas5P())
oplt1 = plot(osol1; title = "No Oscillation")
oplt2 = plot(osol2; title = "Oscillation")

plot(oplt1, oplt2; lw = 3, layout = (2,1), size = (800,600))

The repressilator

The repressilator was introduced in Elowitz & Leibler (2000) as a simple system that can generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in Escherichia col). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) whose production rates are (repressing) Hill functions.

using Catalyst
repressilator = @reaction_network begin
    hillr(Z,v,K,n), ∅ --> X
    hillr(X,v,K,n), ∅ --> Y
    hillr(Y,v,K,n), ∅ --> Z
    d, (X, Y, Z) --> ∅
end

\[ \begin{align*} \varnothing &\xrightarrow{\frac{v K^{n}}{K^{n} + Z^{n}}} \mathrm{X} \\ \varnothing &\xrightarrow{\frac{v K^{n}}{K^{n} + X^{n}}} \mathrm{Y} \\ \varnothing &\xrightarrow{\frac{v K^{n}}{K^{n} + Y^{n}}} \mathrm{Z} \\ \mathrm{X} &\xrightarrow{d} \varnothing \\ \mathrm{Y} &\xrightarrow{d} \varnothing \\ \mathrm{Z} &\xrightarrow{d} \varnothing \end{align*} \]

Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of noise-induced oscillation.

using OrdinaryDiffEq, StochasticDiffEq, Plots
u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
tspan = (0., 200.)
ps1 = [:v => 10.0, :K => 20.0, :n => 3, :d => 0.1]
ps2 = [:v => 10.0, :K => 50.0, :n => 3, :d => 0.1]

oprob1 = ODEProblem(repressilator, u0, tspan, ps1)
oprob2 = ODEProblem(repressilator, u0, tspan, ps2)
osol1  = solve(oprob1, Tsit5())
osol2  = solve(oprob2, Tsit5())
oplt1 = plot(osol1; title = "Oscillation (ODE, K = 20)")
oplt2 = plot(osol2; title = "No oscillation (ODE, K = 50)")

sprob1 = SDEProblem(repressilator, u0, tspan, ps1)
sprob2 = SDEProblem(repressilator, u0, tspan, ps2)
ssol1  = solve(sprob1, ImplicitEM())
ssol2  = solve(sprob2, ImplicitEM())
splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)")
splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")

plot(oplt1, oplt2, splt1, splt2; lw = 2, layout = (2,2), size = (800,600))

The Willamowski–Rössler model

The Willamowski–Rössler model was introduced in Willamowski & Rössler (1979) as an example of a simple CRN model which exhibits chaotic behaviours. This means that small changes in initial conditions can produce relatively large changes in the system's trajectory.

using Catalyst
wr_model = @reaction_network begin
    k1, 2X --> 3X
    k2, X --> 2X
    k3, Z + 2X --> 2Z
    k4, Y + X --> 2Y
    k5, Y --> ∅
    k6, 2Z --> ∅
    k7, Z --> ∅
end

\[ \begin{align*} 2 \mathrm{X} &\xrightleftharpoons[k2]{k1} 3 \mathrm{X} \\ \mathrm{Z} + 2 \mathrm{X} &\xrightarrow{k3} 2 \mathrm{Z} \\ \mathrm{Y} + \mathrm{X} &\xrightarrow{k4} 2 \mathrm{Y} \\ \mathrm{Y} &\xrightarrow{k5} \varnothing \\ 2 \mathrm{Z} &\xrightarrow{k6} \varnothing \\ \mathrm{Z} &\xrightarrow{k7} \varnothing \end{align*} \]

Here we first simulate the model for a single initial condition, showing in both time-state space and phase space how it reaches a strange attractor.

using OrdinaryDiffEq, Plots
u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
tspan = (0.0, 50.0)
p = [:k1 => 2.1, :k2 => 0.7, :k3 => 2.9, :k4 => 1.1, :k5 => 1.0, :k6 => 0.5, :k7 => 2.7]
oprob = ODEProblem(wr_model, u0, tspan, p)
sol = solve(oprob, Rodas5P())

plt1 = plot(sol; title = "Time-state space")
plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")
plot(plt1, plt2; layout = (1,2), size = (800,400))